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DNS and LES of a shear- free mixing layer 

By B. Knaepen, O. Debliquy' and D. Carati f 


1. Introduction 

The shear-free mixing layer represents one of the simplest inhomogeneous flows. It 
consists of two ‘distinct’ homogeneous regions of different turbulent kinetic energy in- 
teracting through a layer of rapid transition. The layer is said to be shear-free since the 
two homogeneous regions have no relative velocity. While flows encountered in nature 
or industrial applications are more often not devoid of shear, the study of the shear- free 
mixing layer is nevertheless useful since it allows the mixing properties of turbulence to 
be examined in a simplified environment. Indeed, turbulent transport properties in shear 
flows are more difficult to track since they can be overwhelmed by production sources 
originating from gradients in the mean velocity. 

The shear- free mixing layer has already received attention in the past, both from the 
experimental and the numerical point of view. The first experimental study of the shear- 
free mixing layer is due to Glibert (1980). The flow was obtained by forcing a stream 
through a grid with two different mesh spacings. The two sides of the grid have however 
equal solidity resulting in an outgoing shear-free flow. In this first experimental study, the 
author mainly concentrated his attention on the downstream evolution of the spreading- 
rate parameter, which is a measure of the thickness of the mixing-layer. A more extensive 
experimental study was later performed by Veeravalli & Warhaft (1987, 1989). Owing to 
a different experimental setup, the authors achieved a higher ratio between the energies 
characterizing the two sides of the flow and this resulted in the observation of large-scale 
intermittency in the flow. The Veeravalli & Warhaft (1989) data will be used here as the 
main benchmark for the present study since it contains a detailed documentation of the 
flow characteristics we will be examining. From the numerical point of view, the shear-free 
mixing layer was studied in Briggs et al. (1996) using DNS. In that article the authors 
also used the data of Veeravalli & Warhaft (1989) as the point of comparison. Their 
simulations were done using a spectral code with a resolution of 128 3 Fourier modes. 
The microscale Reynolds numbers reached in the low- and high-energy homogeneous 
regions were, respectively, 11 and 22.5, which are roughly half of the values reported 
in the Veeravalli & Warhaft (1989) experiment (note that Briggs et al. (1996) use a 
different definition of the microscale Reynolds number than the one in Veeravalli & 
Warhaft (1989) and here). However, when properly nondimensionalized, they were able 
to reproduce satisfactorily the turbulence statistics of the flow. 

The purpose of this work is twofold. First, given the computational resources available 
today, it is possible to reach, using DNS, higher Reynolds numbers than in Briggs et al. 
(1996). In the present study, the microscale Reynolds numbers reached in the low- and 
high-energy homogeneous regions are, respectively, 32 and 69. The results reported earlier 
can thus be complemented and their robustness in the presence of increased turbulence 
studied. The second aim of this work is to perform a detailed and documented LES of 
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Figure 1 . Graphical representation of the shear-free mixing layer 

the shear-free mixing layer. In that respect, the creation of a DNS database at higher 
Reynolds number is necessary in order to make meaningful LES assessments. From the 
point of view of LES, the shear- free mixing-layer is interesting since it allows one to 
test how traditional LES models perform in the presence of an inhomogeneity without 
having to deal with difficult numerical issues. Indeed, as argued in Briggs et al. (1996), 
it is possible to use a spectral code to study the shear-free mixing layer and one can thus 
focus on the accuracy of the modelling while avoiding contamination of the results by 
commutation errors etc. 

This paper is organized as follows. First we detail the initialization procedure used in 
the simulation. Since the flow is not statistically stationary, this initialization procedure 
has a fairly strong influence on the evolution. Although we will focus here on the shear- 
free mixing layer, the method proposed in the present work can easily be used for other 
flows with one inhomogeneous direction. The next section of the article is devoted to 
the description of the DNS. All the relevant parameters are listed and comparison with 
the Veeravalli & Warhaft (1989) experiment is performed. The section on the LES of the 
shear-free mixing layer follows. A detailed comparison between the filtered DNS data and 
the LES predictions is presented. It is shown that simple eddy viscosity models perform 
very well for the present test case, most probably because the flow seems to be almost 
isotropic in the small-scale range that is not resolved by the LES. 


2. Initialization of the flow 

From the numerical point of view, one of the most appealing properties of the shear- 
free mixing layer is the possibility of simulating this flow with a purely spectral three 
dimensional code. Indeed, periodicity can be enforced by considering a second mixing- 
layer, which performs the ‘reverse’ transition compared to the first one. The situation is 
depicted in Fig. 1. This also has the advantage that results gathered from the two mixing 
layers can be averaged to improve the statistics. This possibility will be systematically 
exploited in the presentation of numerical results in the following sections. 

Given a 3D spectral code, the non-trivial part of the simulation is then to build a 
suitable initial condition that mimics the mixing layer. Indeed, because the decaying 
mixing layer is not statistically stationary, it is not acceptable to wait until the initial 
condition is forgotten by the flow due to the stochastic nature of turbulence. Indeed, 
the simulation will remain quite strongly influenced by the initial state of the velocity 
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field. To proceed, it is necessary to introduce a few definitions and notations. The Fourier 
modes associated with the velocity field Ui(x, y , z) will be denoted Ui(k x , k y . k z ) and are 
defined by 

Oi ( k x , ky , kz) — ^ ( Ui (x,y,z)e ■ ( 2 . 1 ) 

X 

Since the mixing layer is homogeneous in two directions it is also convenient to consider 
2D Fourier transforms. Here we take the y direction as the inhomogeneous direction and 
define the 2D Fourier modes Ui(k x ,y,k z ) by 

Ui{k x , y, k z ) = ^2 u i( x > V’ z)e - * kj -' Xj ‘ , (2.2) 

X_L 

where x_l = (x,z) and kj^ = ( k x ,k z ). For convenience, we will also adopt the following 
notations: u x = u, u y = v and u z = w. When the flow is homogeneous and isotropic, a 
common way to initialize the modes Ui(k x ,k y ,k z ) is to fix their amplitudes to match a 
given energy spectra E(k) and assign them random phases in such a way that continuity 
is enforced (Rogallo 1981). One then has, 

{\ui(k x ,k y ,k z )\ 2 ) = A 2 (k), with E(k) = 2Trk 2 A 2 (k) (2.3) 

and k 2 = k 2 + k 2 + k 2 . For the case at hand, we can adopt a similar strategy but consider 
instead the 2D spectra in each plane perpendicular to the direction of inhomogeneity. 
Indeed, in those planes the flow is assumed to be homogeneous and isotropic. We thus 
initialize our flow by imposing the following constraints on the 2D Fourier modes (the 
assignment of the random phases and the treatment of continuity will be described be- 
low), 

< \ui(k x ,y, k z )\ 2 >= B 2 (k±,y), with E(kj_,y) = Trk±B 2 (k±,y) (2.4) 

and k 2 ± = k 2 + k 2 . In the above equation, E(k±,y) is the energy spectra of velocity field 
in (x, z) plane. The arbitrary part remaining is the choice of the function B 2 (k±,y). For 
homogeneous isotropic flows, it is trivial to relate the 2D amplitudes to the 3D amplitudes 
(using Parseval’s theorem): 

B 2 (k ± , y) = j dk v A 2 (k 2 + k\). (2.5) 

Note that as expected, B 2 (k±,y) is independent of y for homogeneous flows. In the 
case of the shear-free mixing layer we will choose an amplitude function A(k) for each 
homogeneous region and compute the corresponding functions B(k±,y) using (2.5). If 
the 2D amplitudes functions in the high energy and low energy regions are respectively 
denoted Bn{k±,y) and BL(k±,y), we then define the complete 2D amplitudes function 
for the shear- free mixing layer as, 

B ML (k±,y) = (1 - f(y))B L (k ± ,y) + f(y)B H (k ±1 y). (2.6) 

The function f{y) is equal to 0 in the low energy region and equal to 1 in the high 
energy region; inside the mixing layers it varies smoothly from 0 to 1. The complete 
initialization procedure is as follows. First we initialize our 3D Fourier modes using the 
procedure of Rogallo (1981). The 3D energy spectra used here is taken from the high- 
energy homogeneous region of the mixing layer. This ensures that the Fourier modes 
Ui(k x , k y , k z ) satisfy the continuity equation. The 3D Fourier modes are then transformed 
to Ui(k Xl y, k z ) using a ID Fourier transform. At this point their 2D amplitudes are 
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measured and the modes are rescaled in order to match the prescribed 2D amplitudes 
given by (2.6): 


Ui(k x ,y,k z ) -> Ui(k x ,y,k z ) = Ui(k x ,y,k z )* 


B ML (k±,y) 

\ui{k x ,y,k z )\ 2 ' 


(2.7) 


The u((fc x; y, k z ) are then transformed back into 3D Fourier modes, u( (fc x , k y , k z ). By 
performing the transformation (2.7), one of course destroys the continuity property of the 
initial field and it has to be recovered by projecting the u((fc x , k y , k z ) onto a divergence- 
free field: 


Ui(k x , ky, k z ^j > U} {k X i ky , k z ) — 


kikj 

6 « ~ V 


U j (k x , ky , kz') . 


( 2 . 8 ) 


This in turn partly undoes the prescription of the 2D amplitudes. Fortunately, by iter- 
ating the transformations (2.7) and (2.8) one converges to a velocity field that has 2D 
amplitudes arbitrarily close to Bml(&_l> 2/) and which satisfies continuity. At this stage, 
the flow still has random phases. In order to correct this problem, we have time evolved 
the flow until the global skewness of the velocity derivatives reached a converged value. 
Between each time step, the mixing layer was rebuilt using (2.7) and (2.8) in order to 
retain the desired 2D amplitudes profile. After that, we stopped the rebuilding procedure 
and let the flow decay freely. 


3. DNS results 


3.1. Parameters of the simulation 

The choice of parameters for the DNS was mainly guided by the following considerations. 
In order to have an experimental reference to compare with, the parameters of the DNS 
have been chosen to match as closely as possible those from the Veeravalli & Warhaft 
(1989) experiment performed using the 3 : 1 perforated plate. Of course, since numerical 
capabilities are not unlimited, some compromises had to be made. The most important 
restriction in the present study is the ability to adequately resolve the high-energy region 
of the flow. Given this constraint, the initial 3D spectra of the homogeneous high-energy 
region En(k) was chosen to match the spectra measured in the Comte-Bellot & Corrsin 
(1971) experiment at stage 1. This experimental spectra was fitted with the following 
function, 


E h {k) 


ak 4 _ hk a 
(k 4 + q A ) 1+a 6 


(3.1) 


which contains several parameters a, q , b, a and (3. This fairly complicated function has 
been chosen because it allows an easy fit of various properties of the energy spectrum. For 
instance, the parameters b and /3 can be used for characterizing the viscous range of the 
spectrum. The parameters q and a determine the energy peak and the transition between 
the energy containing scales and the viscous range. Finally, a determines the total energy. 
The function (3.1) does not allow one to derive analytical expressions for the total energy 
and the total dissipation in terms of the parameters a, q , b, a and (3. It is thus not possible 
to express these parameters in terms of simple global experimental data and we have not 
found a systematic procedure for prescribing them. It has however been observed that 
the following set of parameters, a = 10.6, q = 1.5, b = 0.02, a = 1.233, (3 = 1.1, allows an 
almost perfect fit to the Comte-Bellot & Corrsin spectrum. Of course, the value of some 
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Figure 2. 3D contour plot of the energy density for the initial velocity field. The labels indicated 
on the figure correspond to grid-points. The data field was sampled every four grid-points in 
order to produce the graph. 


of these parameters might depend on the units chosen to perform the simulation. This 
is not really an issue since the time scale and the length scale can be seen as entirely 
defined by the computational domain size l x = 2ir , l y = 47r, l z = 2n and by the viscosity, 
chosen here to have the numerical value of 0.006. In the Veeravalli & Warhaft (1989) 
experiment the ratio of energy between the two homogeneous regions is 6.27 while the 
ratio of dissipation is 7.28. These two ratios can be reproduced well by choosing the 
spectra of the low energy homogeneous region to be of the same form as (3.1) but with 
the following parameters: a = 2.74, q = 3.33,6 = 0.027, a = 1.233, (3 = 1.1. Furthermore, 
with the above choices of parameters, the maxima of the two spectra are separated by a 
ratio that matches the inverse ratio of the initial integral length-scales in the Veeravalli 
& Warhaft (1989) experiment between the low- and high-energy regions. Accordingly, 
the differences in typical sizes of the large-scale structures on both sides of the mixing 
layer are also reproduced. From these definitions it is possible to compute numerically 
the two functions BL(k±,y) and Bn{k±,y) needed in (2.6) using (2.5). The initialization 
procedure is then fully defined if the smoothing function f{y) is prescribed. Here, the 
following choice has been adopted: 


f(y) = < 


0 

l/2(sin(127 r fa ^ y y/4) ) + 1) 

1 

l/2(sin(127r to ~f I,/3) ) + 1) 

0 


if 0 < y < 5/ y /24, 
if 5Z y /24 < y < 7Z y /24, 
if 7l y /24 <y< 17Z y /24, 
if 17ly/2A <y< 19/ y /24, 
if 19Z y /24 <y<l y . 


(3.2) 


With this smoothing function, the high-energy region and the combined low-energy 
regions have the same length. Both are five times larger than each of the two mixing 
layers. As an illustration, the contour plot of the initial energy density is shown in Fig. 2. 
As far as the numerics are concerned, our DNS was performed using a pseudo-spectral 
dealiased code. The grid resolution adopted consists of 512 x 1024 x 512 points. The 
higher resolution dimension being the direction of inhomogeneity, taken here to be y. 


3.2. Kinetic energy diagnostics 

One of the major motivations of this study is to investigate the effect of inhomogeneities 
on the properties of the turbulence. Since the direction of inhomogeneity is along y, it 
is convenient to present statistics obtained by averaging over the x and z directions. 
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Figure 3. (a) Energy and (b) Dissipation rate profiles across the mixing layer calculated from 
(3.3) and (3.4); : t* = 0; : t* = 0.56; : t* = 1.51. 

For instance, the kinetic energy and dissipation rate profiles are calculated from the 
expressions: 

E (y) = ^M x )l 2 = — ^k( x )l 2 ’ ( 3 - 3 ) 

Z 'f'x^z x z ^ 

e{y) = 2 v Sij(x)Sij(x) , (3.4) 

where SV y -(x) = (diUj+djUi)/ 2. Here, and in the rest of this paper, the overbar 7TT denotes 
averaging over the planes perpendicular to the direction of inhomogeneity. Figures 3(a) - 
3(b) represent the profiles of the kinetic energy and the dissipation rate at three different 
times in the simulation. Time has been normalized using the initial eddy turnover time, 
t* = teo/ko where t is the dimensional time, ko is the initial average turbulent kinetic 
energy and eo is the initial average dissipation rate. The Figs. 3 demonstrate that as 
the decay proceeds, the mixing-layer widens but the homogeneous regions remain largely 
discernible. Fig. 4(a) shows the temporal decay of the average energy in the high-energy 
and low-energy homogeneous regions. Assuming an asymptotic power-law decay E(t) ~ 
t n , a decay exponent of n = —1.3 is found in the high-energy homogeneous region while 
in the low-energy region the decay exponent is n = —1.1 (the global energy decay has a 
decay exponent of n = —1.3). These decay rates are compatible with the DNS of Briggs 
et al. (1996) for which a fc 4 low wave-vector energy spectrum was adopted as in the 
present simulation. In Fig. 4(b) the 2D spectra defined in (2.4) are presented to confirm 
that the flow is sufficiently resolved. As was observed in Briggs et al. (1996), the energy 
of the high wave-numbers decays faster with time than the energy of the low wave- 
numbers. It is noted that for the discretization used and the times considered, the energy 
peaks remain at a constant wave-vector in the homogeneous regions. The strong influence 
from the homogeneous layers seems to induce a shift of this energy peak in the mixing 
layer towards higher wave- vectors, though this effect remains moderate for the times 
considered. 


3.3. Intermittency 

The skewness and kurtosis of a velocity component Uj are respectively defined as, 


o _ K 

— 


(«i) ! 


K u = -=^-. 

(«f ) 2 


(3.5) 
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Figure 4. (a) Energy decays of the high ( ) and low ( ) energy homogeneous regions; 

( ) represents the global energy decay. All curves have been normalized using their initial 

value of ko. (b) 2D spectra E(k±,y) (defined by (2.4)) computed inside the high (H) and low 

(L) homogeneous regions and the mixing-layer (M) for three different times; : t* = 0; 

: t* = 0.56; : t* = 1.51. 

For homogeneous isotropic turbulence, measures of these quantities show that they are 
very close to those calculated for a Gaussian signal, i.e. , 5 = 0 and K = 3. 

The skewness profile of v is shown in Fig. 5(a). The y direction has been normalized 
by the half- width li/ 2 of the mixing layer and centered around the inflection point of 
the variance of v. Other skewness and kurtosis diagnostics are normalized in a similar 
fashion (see Veeravalli & Warhaft (1989) for the details of this normalization procedure). 
The skewness profile of v exhibits a sharp deviation from the Gaussian value around the 
location of the mixing layer. As described in Veeravalli & Warhaft (1989), this behaviour 
is attributed to the intermittent penetration into the low-energy region of structures 
originating from the high-energy region (a similar penetration of structures from the 
low-energy region into the high-energy region is certainly also happening but is however 
a lot less frequent). Agreement with experimental data from Veeravalli & Warhaft (1989) 
is very good both in terms of the location of the peak and its amplitude. For symmetry 
reasons, it is expected that the skewness of u and w should remain close to zero. Up 
to statistical deviations, this is confirmed in our simulation (although not illustrated in 
this paper). Kurtosis profiles of the velocity components are shown in Fig. 5(b,c,d) and 
display deviations from the Gaussian value of 3 again around the location of the mixing 
layer. In the Kurtosis of the u, we observe an unexpectedly high peak in the profile 
(compared to the experimental data). At this point we attribute this issue to the fact 
that the profile was computed from a single realization of the flow (although averaging 
in the ( x , z ) planes was performed). We also note, however, an initial small peak in this 
Kurtosis component present in the initial condition (in contrast to the Kurtosis of the 
other two velocity components) but it is not clear whether or not it was amplified during 
the evolution of the flow and by what mechanism. The other two Kurtosis profiles agree 
very well with the experimental data again both in terms of the location of the peaks 
and their amplitudes. Both the numerical simulation and the experiments indicate that 
the deviations from the Gaussian values for S and K occur on the low energy side of the 
inflection points. This supports the idea that these deviations result from the more likely 
penetration of intermittent structures from the high-energy region into the low-energy 
region. Finally, it must be stressed that the initialization procedure presented in section 
2, though quite sophisticated, does not allow imposition of the initial skewness or the 
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Figure 5. (a) Skewness of v; (b) Kurtosis of u; (c) Kurtosis of v; (d) Kurtosis of w; : 

t* = 0; : t* = 0.56; : t* = 1.51; o: experimental data from Veeravalli & Warhaft 

(1989). At t* = 0 the numerical curves are close to their Gaussian value since the penetration 
mechanism described in the text has not occurred yet. 


kurtosis profiles in the region of the mixing layer (Gaussian values are in fact observed 
everywhere for the initial profiles as shown in Fig. 5 except for a small deviation for K u 
as mentioned above). The fact that the DNS later reproduces the experimental profiles 
observed in the mixing layer indicate that the transport mechanisms are successfully 
resolved. 


4. LES of the shear-free mixing-layer 

4.1. Notations and conventions - LES model 
Starting from the Navier-Stokes equations, one obtains the LES equations (4.1) by ap- 
plying a filter, here denoted ~ (since our code is spectral, we will only consider spectral 
cut-offs for the LES filter): 

d t Ui + d jiujUi) = —dip + vAui - djfij. (4.1) 

The unknown subgrid-scale stress (SGS) tensor t ij = UiUj — UiUj needs to be modelled 
in terms of Ui in order to close (4.1). In this work, we will use for r,j a model proposed 
in Wong & Lilly (1994) and further studied in Carati et al. (1995) and in Dantinne et al. 
(1998). This model, which can be considered as a variant of the dynamic Smagorinski 
model (Smagorinsky 1963; Germano et al. 1991; Lilly 1992; Germano 1992), has been 
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shown to perform very well in the context of homogeneous isotropic turbulence and its 
predictions are extremely close to the dynamic Smagorinski model. The advantage of 
this model rests upon the ease with which its dynamic version can be implemented. The 
definition of the model is the following: 

1 . ~ 4 ~ 

X'ij — ft /,: djj — 2C A 3 Sij, (4-2) 

where Sij = \{diUj + djiif) is the resolved strain tensor and A is the LES filter width. 
The dimensional parameter C is evaluated by introducing a second (coarser) filter — 
(the test filter) and using the dynamic procedure: 


C = 


2(A 4 /3 _ A 4 / 3 ) 


Lij Sij > 
< Sij Sij 


(4.3) 


where L,y = — uftij is the Leonard tensor (note that we have systematically used the 

property ~ = — valid for spectral cut-offs). As in the dynamic Smagorinski model, the 
only free parameter available is the ratio of the spectral cut-offs: A/A. In the following 
discussion this ratio will be assumed to be equal to 2. For homogeneous isotropic turbu- 
lence, the averages (• • • ) in (4.3) are obtained by averaging over the whole computational 
domain. The idea is of course that, since turbulence is homogeneous, the constant C 
should be statistically independent of the position. For inhomogeneous flows in one di- 
rection, like the shear-free mixing layer or the channel flow, dependence on the direction 
of inhomogeneity is introduced by averaging quantities only over the other two homoge- 
neous directions. This is justified only if the flow is not too inhomogeneous. The dynamic 
coefficient C then depends explicitly on the inhomogeneous direction: C = C(y). Our 
LES condition was obtained by filtering (with a sharp spectral cut-off) the initial DNS 
field down to 32 x 64 x 32 modes. Thus, only A of the gridpoints are retained in each 
direction for the LES, meaning that there is one LES grid point for about 4 000 DNS 
grid points. The box-size is unchanged and remains 27r x 47t x 27t. In our study we have 
also included a run obtained at LES resolution but with no subgrid-scale stress model to 
emphasize the effect of the model in the LES simulation. For comparison we have filtered 
the DNS fields stored during the simulation down to 32 x 64 x 32 modes. 


4.2. Comparison of the filtered DNS and the LES 

Figure 6(a) represents the temporal evolution of the normalized global kinetic energy 
E/Eq. From the graph it is clear that the simulation with LES model does a very good 
job in reproducing this diagnostic, while the run with no model does not. Likewise, Figs. 
6(b) and 6(c) show that the 2D spectra with LES modelling are in good agreement with 
DNS data. 2D spectra gathered from the run with no model exhibit a strong pile up of 
energy in the high wave-number side of spectra, indicating that the flow is not adequately 
resolved in that case. 

To further illustrate the decay of the kinetic energy, we display in Fig. 7(a) and 7(b) 
the profiles of the kinetic energy at two different times. This is of course a more sensitive 
diagnostic since it retains information about the inhomogeneity of the flow. Both figures 
again reveal a good agreement between the DNS and LES runs and a poor behaviour 
of the run without modelling. Variance profiles of u, v and w have been examined (not 
shown) and again the LES matches the DNS very nicely, whereas the no model simulation 
performs poorly. For our assessment of LES we have retained the same intermittency 
diagnostics described earlier: the skewness and the kurtosis of the velocity components. 




Figure 6. (a) Evolution with time of the normalized global kinetic energy E/Eq\ (b) and (c) 
2D spectra E(k±,y) (defined by (2.4)) computed inside the high (H) and low (L) homogeneous 
regions for time t* = 1.51. In both figures, : filtered DNS; : LES; : no model. 




Figure 7. Energy profile across the mixing-layer calculated from (3.3) at times (a) t* = 0.56, 
(b) t* = 1.51. In both figures, : filtered DNS; : LES; : no model. 


A sample of those quantities is displayed in Fig. 8. Once again, the LES run produces 
results which compare well with the filtered DNS data. Surprising at first, the run without 
subgrid-scale stress model also performs quite well. However, since the intermittency is 
attributed to large-scale structures penetrating the low energy region from the high 
energy region this is to be expected. Indeed, looking at the spectra displayed in Fig. 6(c), 
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Figure 8. (a) Skewness of v ; (b) Kurtosis of u ; (c) Kurtosis of v; (d) Kurtosis of w. All curves 
computed at t* = 1.51. : filtered DNS; : LES; : no model. 


we see that the run without subgrid-scale stress model is not ill-behaved for modes at 
the lowest wave-numbers. 


5. Conclusions 

This study of the shear-free mixing layer reveals that DNS is capable of reproducing 
most of the aspects of the experimental database that have been produced by Veeravalli 
& Warhaft (1989). Energy, dissipation rate and velocity variance profiles are accurately 
reproduced. Also, departure from Gaussianity inside the mixing layer revealed by the 
measurement of the skewness and the kurtosis of the velocity field are also in excellent 
agreement with the experimental observation. This is particularly satisfactory since inside 
the mixing layer the initial values of these quantities cannot be prescribed by the initial- 
ization procedure and are entirely produced at later times by the transfer mechanisms 
simulated by the DNS. It is also observed that, in the absence of shear, no significant 
anisotropy is observed in this flow. This is consistent with previous numerical results and 
strongly supports the use of eddy-viscosity type models for the LES of such flows. As a 
consequence, it is not surprising that comparisons between the predictions of LES using 
such eddy- viscosity models and DNS show very good agreement. 
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